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The process of protein synthesis in biological systems resembles a one dimensional driven lattice 
gas in which the particles (ribosomes) have spatial extent, covering more than one lattice site. 
Realistic, nonuniform gene sequences lead to quenched disorder in the particle hopping rates. We 
study the totally asymmetric exclusion process with large particles and quenched disorder via several 
mean field approaches and compare the mean field results with Monte Carlo simulations. Mean 
field equations obtained from the literature are found to be reasonably effective in describing this 
system. A numerical technique is developed for computing the particle current rapidly. The mean 
field approach is extended to include two-point correlations between adjacent sites. The two-point 
results are found to match Monte Carlo simulations more closely. 



PACS numbers: 82.39.-k, 05.10.-a 
Introduction 



The process of protein synthesis, called "translation," 
] can be modeled using a driven lattice gas in one dimen- 
sion H H H 0. This type of model is well understood 
when all particle hopping rates are uniform, but a model 
for the real biological process requires nonuniform par- 
ticle hopping rates. Direct Monte Carlo simulation of 
such models is possible when only a few genes are in- 
volved. However, it is desirable to perform large-scale 
[ simulations to fit translation models to experimental data 
collected for many genes simultaneously (e.g., data from 
; 13). For this purpose, Monte Carlo approaches would be 
computationally too slow. Therefore, other analytical or 
computational methods are needed. 

In this paper, we address the issue of quenched disor- 
der (site-dependent hopping rates) in a driven lattice gas 

' model for translation. The paper is organized as follows. 
The model is first described and its connection to bio- 
logical protein synthesis explained. Section ITU contains a 
brief summary of known results. In Section lllll we use 
a simple coarse-grained approach to obtain crude, ap- 
proximate solutions. Section llVI treats our central topic: 
application of a mean field method ll, '5| to the problem 

' of quenched disorder. Analytical and computational re- 
sults are presented. In Section O we extend the mean 

' field approach to include two-point correlations for better 
accuracy. We close with a brief summary and discussion 
of how these methods may be applied to problems of in- 
terest, such as fitting translation models to experimental 
data. 
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I. DESCRIPTION OF MODEL 

We focus on translation in prokaryotes, particularly 
Escherichia coli, because of its relative simplicity. The 
process involves the synthesis of specific proteins based 
on the sequence of messenger RNA (mRNA) molecules. 
The mechanism consists of ribosomes "reading" the 
codons of mRNA as the ribosomes move along an mRNA 
chain, and the recruitment and assembly of amino acids 
(appropriate to the codons being read) to form a pro- 
tein. (See, e.g., 0, for more details.) This process is of- 
ten described as three steps: initiation, where ribosomes 
attach themselves, one at a time, at the "start" end of 
the mRNA; elongation, where the ribosomes move down 
the chain in a series of steps; and termination, where 
they detach at the "stop" codon. Since ribosomes can- 
not overlap, their dynamics is subject to the excluded vol- 
ume constraint. Apart from being impeded by another 
ribosome (steric hinderance), a ribosome cannot move 
until the arrival of an appropriate transfer RNA, carry- 
ing the appropriate amino acid (a combination known as 
aminoacyl-tRNA, or aa-tRNA). Thus, the relative abun- 
dances of the approximately 60 types '7| of aa-tRNA have 
significant effects on the elongation rate. Assuming re- 
actant availabilities in a cell are in their steady state, 
with a time-independent concentration of ribosomes and 
aa-tRNA, there would be an approximately steady cur- 
rent of ribosomes moving along the mRNA, resulting in 
a specific production rate of this particular protein. Our 
goal is the prediction of the protein production rates for 
various mRNA's, as a function of the concentration of 
ribosomes and aa-tRNA's. 

The process of translation is well suited to modeling 
using a driven lattice gas in one dimension. Particles 
enter at some rate on one end of a chain of discrete lat- 
tice sites (initiation), then hop down the chain one site 
at a time with another rate or set of rates (elongation), 
and finally exit the chain at the other end (termination). 
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The excluded volume constraint is implemented by in- 
suring the spacing between ribosomes is no less than 12 
sites, the approximate number of codons that a ribosome 
blocks 8,1)]. Quenched disorder arises because there is 
non-uniformity in the hopping (elongation) rates along 
the chain. This effect occurs because at each codon, a 
ribosome has to "wait" for the appropriate aa-tRNA be- 
fore continuing, and the various aa-tRNA's are present 
in nonequal abundances. 

The model we employ is the same as in We model 
an mRNA with N codons as a chain of sites, each of 
which is labeled by i. The first and last sites, i = 1,N, 
are associated with the start and stop codons, respec- 
tively. At any time, attached to the mRNA are M ri- 
bosomes (particles). Being a large complex of molecules, 
each ribosome will cover £ sites (codons), with £ = 12 
typically [1I3. Any site may be covered by a single ri- 
bosome or none. In case of the latter, we will refer to the 
site as "empty" or "occupied by a hole." To locate the 
ribosome, we arbitrarily choose the lowest site covered. 
For example, if the first £ sites are empty, a ribosome can 
bind in an initiation step, and then it is said to be "on 
site i = 1." We define Ui to be the ribosome density at 
site i, where only the left end of the ribosome is counted. 
(In [l|,l3i particles were located by their right end, but 
either choice leads to the same rules of motion.) We also 

define the coverage density pi = Yfs=o ''^i-sj which is the 
probability that site i is covered by some portion of a 
ribosome. 

Next, we specify the dynamics of our model. An at- 
tached ribosome located at site i will move to the next 
site (i + 1) with a rate fc^, provided site i + £ is empty. 
For Monte Carlo simulations, it is convenient to update 
configurations in discrete time units At. In implementing 
the simulations, it is better to use probabilities = fc^Ai, 
so that a ribosome on site i will be moved (or not) with 
probability Pi (or l—pi). We purposefully associate these 
hopping probabilities with a site because a site is asso- 
ciated with a particular codon. Thus, the hopping rate 
from that site may depend on the relative abundance of 
the appropriate aa-tRNA. Apart from these probabili- 
ties, another aspect of our simulations is random sequen- 
tial updating: i.e., during each Monte Carlo step (MCS), 
M + 1 particles are chosen at random, in sequence, to at- 
tempt moves. They are selected from a pool that includes 
the M particles on the lattice plus another unbound par- 
ticle that can initiate if there are £ holes at the beginning 
of the chain. 

In our computational studies, systems begin empty and 
are run for long enough to reach steady state. After 
steady state is attained, data including the current and 
density distribution can be collected. Density data is typ- 
ically collected every 100 MCS. We often use continuous 
time Monte Carlo |l3 because it runs far more quickly 
than and provides the same results as standard Monte 
Carlo. Using continuous time Monte Carlo also avoids 
the need to specify a fixed time step At. 

In our studies of real mRNA sequences, we use gene se- 



quences from Escherichia coli strain MG1655, obtained 
from Elongation rates at each codon are estimated 
using commonly accepted values for the availability of 
tRNA in E. coli T^l . The rate at each codon is assumed 
proportional (with an arbitrary proportionality constant) 
to the availability of the appropriate tRNA that can de- 
code the codon, as in '13]. Corresponding data are not 
available for estimating initiation and termination rates, 
so we choose various rates of interest to study. 



II. SUMMARY OF PREVIOUS RESULTS 

Extensive investigations of the simple totally asymmet- 
ric exclusion process (TASEP, defined as point particles 
hopping with unit rate along a line) with open boundaries 
can be found in the literature. We first consider studies 
of uniform systems. Exact analytic results for the £ — 1 
steady state exist ^M- Systems with extended ob- 
jects {£ > 1) have been less frequently investigated but 
have also been understood from various approaches. Us- 
ing a mean field approach, MacDonald et al. derived 
mean field equations for the average site occupation (n^) 
and considered both closed 1] and open Q systems. In 
the former case, exact solutions were found, leading to 
a current vs. density relation. For the latter, the au- 
thors resorted to numerical solutions to find the phase 
diagram for a variety of initiation and termination rates. 
Lakatos and Chou j3| considered uniform open systems 
with extended objects. Using a discrete Tonks gas par- 
tition function, they derived the current vs. density re- 
lation first presented by MacDonald et al. 0]. Via a 
refined mean field theory, they extended this result to 
predict currents and bulk densities for the open system, 
which they confirmed by Monte Carlo simulations. Fi- 
nally, Shaw et al. 0| used an extremal principle 
based on domain wall theory ^3 to obtain the phase di- 
agram, with currents and bulk densities, for the uniform 
open system. They also found approximate density pro- 
files (related to the coverage density p) from a continuum 
limit. From all of these studies, the £ > 1 phase diagram 
is well known. Depending on the initiation (or injec- 
tion) and termination (or depletion) rates, the system 
will settle into one of three phases. From their dominant 
characteristics, the three phases are known as low den- 
sity, high density, and maximal current. The initiation 
and termination probabilities are often referred to as a 
and (3 in the literature. A phase diagram in this a-/3 
plane has been determined, showing second order transi- 
tions between the maximal current phase and the others, 
as well as a first order transition between the high- and 
low-density regions. 

When disorder is introduced, i.e., not all the pi's are 
equal, then methods for exact analytic approaches fail 
(except in the extremely dilute limit, where only the mo- 
tion of a single particle is of concern ^3)- Indeed, even 
a single slow rate in a closed system poses serious diffi- 
culties 0,113,13. However, Kolomeisky j22| obtained 
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approximate steady state solutions and phase diagrams 
for an open system with a single nonuniform rate in the 
bulk by splitting the system into two smaller systems 
connected by the nonuniform rate. Tripathy and Barma 
p^ l considered a closed system, but with a finite fraction 
of identical slow sites. Based on a combination of Monte 
Carlo simulations and numerical solutions of mean field 
equations, they found current-density relations. Further 
references on TASEP with disorder may be found in a 
recent review ■ All of these studies are restricted to 
£ — 1. Studies of disorder in systems with £ > 1 have 
been fairly limited. Shaw et al. [j| found upper and lower 
bounds for the current in systems with arbitrary disor- 
der. In another work, Shaw et al. [2^ considered an open 
system with i > 1 and a single nonuniform rate in the 
bulk. As was done for £ = 1 (22i] , the system was divided 
into two smaller systems connected by the nonuniform 
rate. Steady state currents and bulk densities to either 
side of the nonuniform site were obtained. 



III. SIMPLE COARSE-GRAINING APPROACH 

We consider briefly an approximate method motivated 
by the method of Shaw et al. . Their particle hopping 
rate in the bulk was 1, except for the nonuniform rate 
q at special site k. Important in their analysis is the 
parameter 



leff 



l/q+{£~l)- 



This parameter appears in the current passing from the 
left sublattice into the right sublattice: 



J = leff 



Pleft (1 ~ Pright) 

£ (1 

Pright H~ Pright 

l£y 



(1) 



where pieft and Pright are the bulk densities in the left 
and right sublattices. We note that q^jf is the same as 
£Kg^k in the notation of where 



i+l-l ^ 
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FIG. 1: Coverage density profile for the ompA gene of E. coli 
when elongation rates are limiting. Bold curve is Monte Carlo 
simulation result, and lighter curve is coarse-grained predic- 
tion from Eq. |3] using the Monte Carlo current. The positive 
root was used to the left of the current-limiting region and 
the negative root to the right. The real part is plotted where 
the predicted p is complex. Elongation rates at each codon 
were assumed proportional to availabilities of corresponding 
tRNA [3. 



It is reasonable to use the positive (high density) root to 
the left of the current-limiting region (where the mini- 
mum Kg occurs) and the negative (low density) root to 
the right. Results of this approach are shown in Fig. 
n for the ompA gene when elongation rates are limit- 
ing (i.e., initiation and termination rates are sufficiently 
large). The value for current J in Eq. |21is taken from 
Monte Carlo simulations of ompA. The agreement be- 
tween the coarse-grained result and Monte Carlo simu- 
lations is good in low density regions but is poorer in 
high density regions because long range correlations be- 
come more important, an effect not captured by coarse- 
graining over the relatively short distance i. Finally, we 
note that Eq. |3 can be used only when the current J 
is known, either from Monte Carlo simulations or from 
some analytical means. 



is a coarse-grained rate for translating £ sites. 
The form of Eq. 2] motivates us to suggest that 



J = £K. 



1 



£ 1 



■P^/i 



(2) 



in regimes in which the coverage density p is slowly vary- 
ing in i. Because p and Kg are both coarse-grained over 
a distance £, a relationship between them is unsurprising. 
Eq. 12 can be solved for pf. 



Pi 



Ke,, + J-J/£ 



± 



(3) 



IV. MEAN FIELD APPROACH 

We next turn to a mean field approach, using equations 
first developed by MacDonald et al. QjQ- The equations 
were applied only to uniform systems, but we will find 
that they can be successfully applied to nonuniform sys- 
tems. Here the location of a particle is determined by the 
location of its left end. The particle density at site i is 
TT-i, and the hole density at site i is 1 — X]f=o '^i-s- For a 
particle to move from site j to z + 1, producing a current, 
site i + £ must be empty. The method is "mean field" 
in the sense that some correlations have been neglected. 
The conditional probability that site i + £ is empty given 
that site i contains a particle is replaced by the condi- 
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tional probability that site i + £ is empty given that site 
i contains cither a particle or a hole. That is, 

^ ' ' P(.-o) + P(.-.) 

« P(?-o I?-?) 

_ P(« - o) + P(o - o) 

~ P(« - o) + P(« - •) + P(o - o) + P(o - •) 

I 



dni 

~dr 

drii 
~dt 

di 

drii 
~dt 

where we use to denote the initiation rate. 

For the steady state solution, time derivatives are set 
to and the current J is introduced. The resulting set 
of equations, 

J = fco(l ~ X!"'') ^^^^ 

s=l 

_ J I - rij+s + rij+t 

for i = 1, . . . , TV - ^ (5b) 
fii = ^ tor N ~ £+1,.. .,N, (5c) 

can be solved iteratively for un to ni if J is specified. 
Then Eq. I5al becomes an initial condition to check for 
consistency to determine whether J has been chosen cor- 
rectly. If Eq. Eal is not satisfied, then J should be ad- 
justed appropriately and the process repeated. 

We present an argument for the validity of this it- 



where we use • to denote a filled site and o to denote an 
empty site. 



The mean field assumption for the conditional proba- 
bility leads to the following equations for the time evolu- 
tion of the {ui}: 



erative approach. First, a physically meaningful solu- 
tion will have particle density rii € (0, 1) and coverage 

density Es=i ^ (0)1) all i. (Endpoints of the 
interval are excluded if there is to be nonzero current 
flow.) Suppose that for some i, J > ki{l ~ Es=i i^i+s)- 
Then from Eq. I5bl > 1 — Ef=i'^i+S' meaning that 
Es=o ^i+s > Ij which is a contradiction. So we see that 
J < ki{l — Ef=i iT-i+s) for all i. (Note that although we 
have not dealt separately with i — N — £ + 1, . . . , N, Eq. 
I5cl is consistent with the previous statement.) Thus, J 
cannot be too large if physical solutions are to be ob- 
tained. 

Next, we show that the densities {rii}, while within 
physical ranges, are increasing functions of J. Consider 
two different J values: Jo with its associated densities 
{ui} and Ji with its densities {m^}, which we assume to 
be in physical ranges. Suppose that Ji > Jq. Clearly 
nii > Tii for « = N — £+1, . . . ,N. Using induction on the 
remaining i, it can be shown from Eq. 15 bl that 



+ s 



Us - kini- ^ 



1 -J2l=int-i+s 



1 - Es=l "■i-l+s + ni-i+i 1 - Es=l ^i+s + ni+i 



ioii^2,...,N~£ 



1 1 - Es=l "Af-^+s 7 
KN-inN-l —I «;Ar-^+l"'W-£+l 

1 - Z/s=l ^N-i+s + riN 

ki-iUi-i — kiUi for i — N~£ + 2,...,N, 

I 



(4) 



rrii - Ui > 



^oi-e:=i' 



rui+s + m 
Jo 



.+t Jo 1 - Es=l 



^« 1 - eLi 



> 0. 
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Therefore, the densities {rii} increase with increasing J. 

Finally, we again consider Eq. |23 The left side in- 
creases monotonically with increasing J, and the right 
decreases monotonically with increasing J, while den- 
sities are in physical ranges. If we follow the itera- 
tive approach, J values that lead to J2l=i "^i+s > 1 
or J > ko (^l — X]s=i arc too large and should be 
decreased. On the other hand, J values that lead to 
J < A:o ^1 — too small and should be in- 

creased. One can start with upper and lower bounds for 
the current (from 0) and use bisection to converge to the 
correct current. If a physical solution exists, it is unique 
and should be found by this method. 

Note that the upper bound for J from Q, 



J < 



E — 



(6) 



for all i, applies also to the mean field equations. This 
can be shown from 



E 



E^ 



H+s 



H+s+t 



TABLE L Currents for real gene sequences of E. coli from 
Monte Carlo (MC) simulations and mean field (MP) calcula- 
tions. Both the original mean field and the two-point mean 
field are included. Units for the currents are arbitrary. Per- 
cent errors for the mean field currents, as compared with the 
Monte Carlo currents, are given in parentheses. Elongation 
rates were assumed proportional to the availability of the ap- 
propriate tRNA for each codon fl^ . Several of the genes 
were chosen to be initiation-limited by making the initiation 
rate 0.78 of the slowest elongation rate. Others were made 
termination-limited by making the termination rate 0.34 (for 
asnS and envY) or 0.52 (for fabG) of the slowest elongation 
rate. The remainder were elongation-limited. Errors in the 
Monte Carlo results are less than 0.001. 

gene limit MC J orig. MP J (%err) 2-pt. MP J (%err) 



adk 


elong. 


0.139 


0.133 (4.3) 


0.137 (1.4) 


cysK 


elong. 


0.120 


0.122 (1.7) 


0.118 (1.7) 


gapA 


elong. 


0.194 


0.191 (1.5) 


0.191 (1.5) 


gluH 


elong. 


0.156 


0.154 (1.3) 


0.158 (1.3) 


aceP 


init. 


0.170 


0.164 (3.5) 


0.166 (2.4) 


err 


init. 


0.172 


0.177 (2.9) 


0.172 (0.0) 


fabD 


init. 


0.114 


0.112 (1.8) 


0.112 (1.8) 


asnS 


term. 


0.114 


0.114 (0.0) 


0.114 (0.0) 


envY 


term. 


0.092 


0.091 (1.1) 


0.091 (1.1) 


fabC 


term. 


0.112 


0.113 (0.9) 


0.112 (0.0) 



< 



s=l 



< 1. 



In practice, computing iterative solutions for {ui} and 
adjusting J by bisection is effective in finding J and 
finding n; values that are low density (downstream of 
the current- limiting region) . Table ^ shows that there 
is fairly good agreement (within 5%) between the mean 
field current and the actual current (from Monte Carlo 
simulations) for various real gene sequences of E. coli. 
However, numerical problems arise in finding rii values 
that are high density (upstream of the current-limiting 
region). For high density solutions, we have observed 
that there generally exists a very narrow range for J, 
with a width less than machine precision, for which the 
smaller J values will fail to satisfy Eq. [Sa| because the 
densities are too small, and for which the larger J values 
will lead to nonphysical densities. An example of this 
phenomenon is shown in Fig. [21 

We next present an argument for why high density 
solutions are associated with numerical difhculties. For 
convenience, we assume that the ki are uniformly 1, and 
we seek uniform density solutions. Eq. I5bl gives an itera- 
tive map for n^. We assume that a fixed point n* exists. 
It will then satisfy 



= J 



l-(^-l)^ 
1 - in* 




PIG. 2: Particle density profiles calculated by iteration of 
mean field equations (Eq. |5J| for a uniform system with ini- 
tiation rate 1, elongation rates 1, and termination rate 0.1. 
The system had A'^ = 150 and ^ = 8. The dark curve is the 
result for current J slightly too large, and the light/diffuse 
curve is the result for J slightly too small. The two curves 
are superimposed for i > 50. The difference between the two 
J values was 2 x 10"^® J. 



We find high density and low density fixed points. 



1 

2e 



i + j{i-i)±^[i + j{i- i)Y - uj 



Suppose the densities n^+i, . . . , rii+i are slightly per- 
turbed from the high density fixed point, so that Uj = 



6 



n+ + 5, where \5\ ^ 1, for j = i + 1, i + £. Then 



m = J- 



l-£n+- 15 



4J 



-l + {£-\)J+ ^J[l + J{£-l)f-AU 



-0 (6^) 



To determine stabiUty of the high density fixed point, 
we consider the 5 prefactor: 

4J 



-1 + - 1) J + ^^[1 + J(^-1)]^-4^J 



For currents in the expected range, 

(cf. Q), a > 1 so that the high density fixed point is 
unstable. A similar argument shows the stability of the 
low density fixed point. Although these ideas are proven 
here only for uniform density solutions, our numerical 
results (such as those in Fig. \^ lead us to believe that 
the nonuniform density case is similar. It appears that 
small numerical imprecisions prevent us from accurately 
finding high density solutions, while low density solutions 
can be easily found. 

Finding steady state high density mean field solutions 
is a nontrivial problem. We have attempted multidimen- 
sional Newton's method approaches to solve Eq. |S1 for 
{ni\ and J, but these have their own difficulties. Con- 
vergence often fails unless the initial guess is very near 
the solution. The most reliable method is to start with 
an empty lattice and integrate Eq. ^ to steady state. 
Although integration is computationally slow, it consis- 
tently produces density profiles that are reasonable and 
similar to the Monte Carlo simulation density profiles. 
Agreement is best in the low density regime, when the 
correlations neglected by the mean field theory are less 
important (data not shown). In the high density regime, 
the mean field results underestimate the correlations be- 
tween adjacent particles. An example of this discrepancy 
is shown for a uniform system in Fig. |3| 



V. MEAN FIELD WITH TWO-POINT 
CORRELATIONS 

To obtain more accurate solutions for density profiles, 
especially in the high density regime, we extend the mean 
field theory to include two-point correlations between ad- 
jacent sites. Variables in the two-point mean field theory 
are densities of bonds, indexed by «, where bond i con- 
nects site i to site « -f 1. There are four types of bond 
densities, which we denote as follows: no,i, hole- hole pairs 
(o— o); n_,_i, particle-hole pairs (•— o); ri^.i, holc-particle 
pairs (o — •); and nx,i, particle-particle pairs (• — •). Ge- 
ometry requires that 
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FIG. 3; Particle density profile for a uniform system with 
initiation rate 1, elongation rates 1, and termination rate 0.1. 
The system had N = 150 and I = 12. Density peaks (every I 
sites) are displayed as symbols: filled squares are the Monte 
Carlo simulation result, open triangles are the prediction from 
the original mean field theory, and open diamonds are the 
prediction from the two-point mean field theory. Non-peak 
densities are displayed as curves; non-peak results from all 
three methods overlap. 



and 

A third equation. 



1 



n_ 



("x,i+s + f^«-,^-^-s), 



(8) 



can be written because each lattice site is occupied by 
either a hole or some part of a particle. However, this 
third equation is linearly dependent on Eqs. [3and|Hl 

We can use Eqs. [3and|Slto eliminate n^^i and Uo^i 
from the problem, so we will write differential equations 
for the time evolution of only the two remaining types 
of densities, n^.i and rix.j- Fluxes into and out of each 
state take the form of a product of the appropriate rate 
constant, the density of the state that may evolve, and 
the conditional probability that adjacent particles and 
holes are arranged appropriately for the evolution to oc- 
cur. For example, a hole-hole pair at bond i will evolve 
to a particle-hole pair at bond i with rate 



I? - o - 



0- 



We make mean field assumptions for the conditional 
probability, similar to that of MacDonald et al. |1| . For 
example. 



Pi 



• — o — o 



(7) 



? - o - o) 

P(« - 0-? I? - 0-?) 

n^.i-i + rioA-i 
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Thus the flux of hole-hole pairs at site i to particle-hole The resulting differential equations for time evolution 
pairs at site i is of the bond densities in the bulk are 

n^.i-i + no ,-1 



— -— = ki+in^i ■ Vki-iUoA '- kiU^i (9a) 

at n^^i+i + riy^^i+i n^_i_i -|- no,i_i 

dn^i .i-i , n— 

— = ki^in^ i Ki^iTix i ■ ■ (yb) 

at ' n^_i_i -f 7io,i-i n^^i+e + nx,i+£ 

We also have boundary conditions, and for convenience, we assume that particles enter the lattice one site at a time, 
so that only the first site must be free for initiation to occur. We also assume that a particle whose right edge is on 
site TV can leave the lattice, freeing the final £ sites. Particle-particle bonds thus cannot exist in the final £ sites. Then 
the boundary conditions are, for initiation, 

dn^,i n^^i+e , , , 

— -r— = /ci+^nx.i ; h Kono.i - fcin^ 1 

at n^s+i + n^.i+i 

dn-^i n^,i+£ 



dt ' ' n^^i+i + nx,i+e 

and, for i = iV - ^ + 1, . . . ,7V, 

at n^^i-i + Uo^i-i 

driY i 



dt 



It is possible to obtain an iterative steady state solution 
for the bond densities, from i = iV to i = 1. However, 
this solution appears to exhibit numerical instabilities in 
the high density regime that are similar to those observed 
with the original mean field theory. We would prefer to 
have a simple method, like the iteration and bisection 
method, for computing the current despite the numerical 
difficulties in computing the densities. However, such a 
method is not apparent. Instead, we begin with an empty 
lattice and directly integrate the differential equations for 
the bond densities (Eq. O until steady state is attained. 
The two-point mean field approach produces both den- 
sities and currents that agree more closely with Monte 
Carlo simulations than did the original (one-point) mean 
field theory. Table |I| compares two-point mean field cur- 
rents with currents from Monte Carlo and the one-point 
mean field theory for real gene sequences. In each case, 
the two-point mean field does as well as or better than 
the original mean field at matching the Monte Carlo cur- 
rents. Fig. compares the two-point mean field density 
profile with that obtained by the other methods, showing 
that the two-point mean field theory successfully incor- 
porates more of the long-range correlations than does the 
one-point mean field theory. 



VI. CONCLUSIONS 

We have considered one-dimensional driven lattice gas 
models with large particles and quenched disorder. Mean 
field theories were found be effective in computing quan- 
tities of interest. The mean field equations originally 
proposed by MacDonald et al. 0, for uniform sys- 
tems were found to work equally well for nonuniform 
systems. An iterative method allowed easy and rapid 
computation of the steady state current through the sys- 
tem. Steady state particle densities were also computed 
by this method, although only when the system was in 
a low-density phase. We have gained some insight into 
the numerical difficulties that arise in obtaining high den- 
sity solutions. Direct integration of differential equations 
for the time evolution of particle densities can always be 
used to find the steady state densities. We found good 
agreement between the mean field current and the Monte 
Carlo current. Agreement between the densities was also 
adequate, though not as good in the high density regime. 

We extended the mean field approach to two-point 
correlations, using similar mean field approximations for 
conditional probabilities. Currents and particle densities 
were obtained more accurately from the two-point mean 
field theory than from the original. 

Although the two-point mean field theory is an im- 
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provement on the original theory, Fig. shows that it 
still does not capture all of the correlations necessary to 
reproduce high density profiles. The theory could be fur- 
ther extended to include three-point correlations (such as 
the density of particle-particle-hole triplets). However, 
the number of independent variables and equations, as 
well as the complexity of the equations, would increase as 
more correlations are added. Also, numerical instabilities 
might still exist so that solutions could be found only by 
integration. We therefore feel that it is not convenient to 
extend the method to include higher-order correlations. 

We conclude that mean field approaches can be effec- 
tive in studying disordered systems. If one is primarily 
interested in the current through the system, this quan- 
tity can be computed rapidly using iteration and bisec- 



tion. We expect the iteration and bisection method to 
be quite valuable in future studies, because the calcu- 
lated protein production rates could be compared to ex- 
perimental data (e.g., data in [^) and used in fitting of 
unknown rate constants. 
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